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ABSTRACT 

A formalism is presented that allows cosmological experiments to be 
tested for consistency, and allows a simple frequentist interpretation of 
the resulting significance levels. As an example of an application, this 
formalism is used to place constraints on bulk flows of galaxies using 
the results of the microwave background anisotropy experiments COBE 
and SP91, and a few simplifying approximations about the experimen- 
tal window functions. It is found that if taken at face value, with the 
quoted errors, the recent detection by Lauer and Postman of a bulk 
flow of 689 km/s on scales of 150/i~^Mpc is inconsistent with SP91 at 
a 95% confidence level within the framework of a Cold Dark Matter 
(CDM) model. The same consistency test is also used to place con- 
straints that are completely model- independent, in the sense that they 
hold for any power spectrum whatsoever — the only assumption being 
that the random fields are Gaussian. It is shown that the resulting 
infinite-dimensional optimization problem reduces to a set of coupled 
non-linear equations that can readily be solved numerically. Applying 
this technique to the above-mentioned example, we find that the Lauer 
and Postman result is inconsistent with SP91 even if no assumptions 
whatsoever are made about the power spectrum. 
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1. INTRODUCTION 

Together with the classical cosmological parameters h, fl, etc., the power spectrum P{k) of cos- 
mological density fluctuations is one of the most sought-for quantities in modern cosmology, vital for 
understanding both the formation of large-scale structure and the fluctuations in the cosmic microwave 
backgroimd radiation (CMB). 

The traditional approach has been to assume some functional form for P{k) (like that predicted 
by the cold dark matter (CDM) scenario, for instance), and then investigate whether the predictions 
of the model are consistent with experimental data or not. The large amounts of data currently 
being produced by new CMB experiments and galaxy surveys, all probing different parts of the power 
spectrum, allow a new and more attractive approach. We can now begin to probe exact shape of 
the function P{k), without making any prior assumptions about P{k). More specifically, we measure 
different weighted averages of the function, the weights being the experimental window functions. 

This new approach is quite timely (Juskiewicz 1993), as there are now many indications that 
the primordial power spectrum may have been more complicated than an n = 1 power law. There 
are several sources of concern about the standard CDM cosmology, with inflation leading to f2 ~ 1 
and a primordial n w 1 Harrison- Zel'dovich power spectrum. Compared to COBE-normalized CDM, 
observational data shows unexpected large-scale bulk flows (Lauer & Postman 1993), too weak density 
correlations on small scales (Maddox et al. 1990), a rather quiet local velocity field (Schlegel et 
al. 1993) and a deficit of hot x-ray clusters (Oukbir & Blanchard 1992). The combined data from 
the COBE DMR (Smoot eA, al. 1992) and the Tenerife anisotropy experiment (Watson et al. 1992) 
point to a spectral index exceeding unity (Watson & Gutierrez de la Cruz 1993) which, if correct, 
cannot be explained by any of the standard inflationary models. The recent possible detections of 
halo gravitational microlensing events (Alcock et al. 1993) give increased credibility to the possibility 
that the dark matter in our galactic halo may be baryonic. If this is indeed the case, models with 
O < 1 and nothing but baryonic dark matter (BDM) (Peebles 1987, Gnedin & Ostriker 1992, Cen, 
Ostriker & Peebles 1993) become rather appealing. However, in contrast to CDM with inflation, BDM 
models do not include a physical mechanism that makes a unique prediction for what the primeval 
power spectrum should be. Rather, the commonly assumed P{k) oc is chosen ad hoc to fit 

observational data. Moreover, for fluctuations near the curvature scale in open universes, where the 
f2 = 1 Fourier modes are replaced by hyperspherical Bessel functions with the curvature radius as a 
built-in length scale, the whole notion of scale-invariance loses its meaning (Kamionowski & Spergel 
1993). 

In summary, it may be advisable to avoid theoretical prejudice as to the shape of the primordial 
power spectrum. In this spirit, we will develop a consistency test that requires no such assumptions 
whatsoever about the form of the power spectrum. This approach was pioneered by Juszkiewitz, 
Gorski and Silk (1987), who developed a formalism for comparing two experiments in a power-spectrum 
independent manner. We generalize this method to the case of more than two experiments, and then 
use the formalism to assess the consistency of three recent observational results: the CMB anisotropy 
measurements made by the COBE Differential Microwave Radiometer (Smoot et al. 1992), the South 
Pole anisotropy experiment (SP91, Gaier et al. 1992), and the measurement of bulk velocity of Abell 
clusters in a 150 Mpc sphere (Lauer and Postman 1993, hereafter LP). 

In Section 2, we develop a formalism for testing cosmological models for consistency. In Section 
3, we apply this formalism to the special case of cold dark matter (CDM) and the LP, SP91 and 
COBE experiments. In Section 4, we solve the variational problem that arises in consistency tests of 
models where we allow arbitrary power-spectra, and apply these results to the LP, SP91 and COBE 
experiments. Section 5 contains a discussion of our results. Finally, the relevant window functions are 
derived in the Appendix. 
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2. CONSISTENCY TESTS FOR COSMOLOGICAL MODELS 

In cosmology, a field where error bars tend to be large, conclusions can depend crucially on 
the probabilistic interpretation of confidence limits. Confusion has sometimes arisen from the fact 
that large-scale measurements of microwave background anisotropies and bulk flows are fraught with 
two quite distinct sources of statistical uncertainty, usually termed experimental noise and cosmic 
variance. In this section, we present a detailed prescription for testing any model for consistency with 
experiments, and discuss the appropriate probabilistic interpretation of this test. By model we will 
mean not merely a model for the underlying physics, which predicts the physical quantities that we 
wish to measure, but also a model for the various experiments. Such a model is allowed to contain 
any number of free parameters. In subsequent sections, we give examples of both a very narrow class 
of models (standard CDM where the only free parameter is the overall normalization of the power 
spectrum), and a wider class of models (gravitational instability with Gaussian adiabatic fluctuations 
in a flat universe with the standard recombination history, the power spectrum being an arbitrary 
function). 

Suppose that we are interested in N physical quantities Ci,...,cn, and have A'' experiments 
El,..., En devised such that the experiment Ei measures the quantity Ci. Let .s; denote the num- 
ber actually obtained by the experiment Ei. Because of experimental noise, cosmic variance, etc., we 
do not expect Sj to exactly equal Cj. Rather, Si is a random variable that will yield different values each 
time the experiment is repeated. By repeating the experiment M times on this planet and averaging 
the results, the uncertainty due to experimental noise can be reduced by a factor ^/M. However, if the 
same experiment were carried out in a number of different horizon volumes throughout the universe 
(or, if we have ergodicity, in an ensemble of universes with different realizations of the underlying 
random field), the results would also be expected to differ. This second source of uncertainty is known 
as cosmic variance. We will treat both of these uncertainties together by simply requiring the model 
to specify the probability distribution for the random variables si. 

Let us assume that the random variables Sj are all independent, so that the joint probability 
distribution is simply the product of the individual probability distributions, which we will denote 
fi{s). This is an excellent approximation for the microwave background and bulk flow experiments 
we will consider. Finally, let Si, ...,sjv denote the numbers actually obtained in one realization of the 
experiments. 

The general procedure for statistical testing will be as follows: 

•First, define a parameter rj that is some sort of measure of how well the observed data Si agree 
with the probability distributions fi, with higher rj corresponding to a better fit. 

•Then compute the probability distribution fniv) of this parameter, either analj^tically or by 
employing Monte-Carlo techniques. 

•Compute the observed value of r], which we will denote fj. 

•Finally, compute the probability P{r] < fj), i.e. the probability of getting as bad agreement as 

we do or worse. 
We will now discuss these four steps in more detail. 

2.1. Choosing a goodness-of-flt parameter 

Obviously, the ability of to reject models at a high level of significance depends crucially on making 
a good choice of goodness-of-fit parameter r/. In the literature, a common choice is the likelihood 
product, i.e. 

N 

Tll^WUiSi). (1) 

i=l 

In this paper, we will instead use the probability product, i.e. the product of the probabilities Pi that 
each of the experiments yield results at least as extreme as observed. Thus if the observed Si is smaller 
than the median of the distribution fi, we have Pi = 2P{si < s,), whereas Si larger than the median 
would give Pi — 2P(sj > Si). The factor of two is present because we want a two-sided test. Thus 
Pj = 1 if Si equals the median, Pi = 2% if Si is at the high 99th percentile, etc. 
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2.2. Its probability distribution 



Apart from the simple interpretation of the probability product rj, it has the advantage that its 
probability distribution can be calculated analytically, and is completely independent of the physics 
of the model — in fact, it depends only on N. We will now give the exact distributions. 

By construction, < ?7 < 1. For N = 1, rj will simply have a uniform distribution, i.e. 

Thus in the general case, rj will be a product of A'' independent uniformly distributed random variables. 
The calculation of the probability distribution for r] is straighforward, and can be found in a number 
of standard texts. The result is 

,^^/lAr^(-l'^^)'^"' ifO<77<l, 



L otherwise. 



2.3. The consistency probability 

The probability P{r} < rf), the probability of getting as bad agreement as we do or worse, is simply 
the cumulative distribution function F^{ff), and the integral can be carried out analytically for any N: 

F.ifj) ^p{n<f,)= fr,{u)du = fje{f)) J2 ^ . (4) 

•'0 i=0 

where 9 is the Heaviside step function, and Fr,{fj) = 1 for f) > 1. Since the product of N numbers 
between zero and one tends to zero as A' ^ oo, it is no surprise that 

F^{fj)^e{fi)f,e-'^^ = e{fi) (5) 

as N ^ oo, i.e. that /,,(77) S{fi). The function Fjj{ri) is plotted in Figure 1, and the values of r) for 
which Frj{fi) = 0.05, 0.01 and 0.001, respectively, are given in Table I for a few A^-values. For example, 
if three experimental results give a goodness-of-fit parameter fj = 0.0002 for some model, then this 
model is ruled out at a confidence level of 99%. Thus if the model where true and the experiments 
where repeated in very many different horizon volumes of the universe, such a low goodness-of-fit value 
would be obtained less than 1% of the time. 



Table 1. Probability product limits: 



Confidence level 


N = 1 


N = 2 


N = 3 


N = 4 


95% 


0.05 


0.0087 


0.0018 


0.00043 


99% 


0.01 


0.0013 


0.00022 


0.000043 


99.9% 


0.001 


0.000098 


0.000013 


0.0000021 



2.4. Ruling out whole classes of models 

If we wish to use the above formalism to test a whole set of models, then we need to solve 
an optimization problem to find the one model in the set for which the consistency probability is 
maximized. For instance, if the family of models under consideration is standard n = 1, F = 0.5 CDM 
(see Section 3), then the only free parameter is the overall normalization constant A. Thus we can 
write the c;onsistency probability as p{A), and use some numerical method to find the normalization 
At, for which p{A) is maximized. After this, the statistical interpretation is clear: if the experiments 
under consideration are carried out in an ensemble of CDM universes, as extreme results as those 
observed will only be obtained at most a fraction p{A*) of the time, whatever the true normalization 
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constant is. Precisely this case will be treated in the next section. For the slightly wider class of models 
consisting of CDM power spectra with arbitrary A, n and F, the resulting optimization problem would 
be a three-dimensional one, and the maximal consistency probability would necessarily satisfy 

p{A„n,,r,) > p{A,, 1,0.5) = p{A,). (6) 

An even more general class of models is the set of all models where the random fields are Gaussian, 
i.e. allowing completely arbitrary power spectra P. In section 4, we will show that the resulting 
infinite-dimensional optimization problem can in be reduced to a succession of two finite-dimensional 
ones. 

3. COLD DARK MATTER CONFRONTS SP91, COBE AND LAUER-POSTMAN 

As an example of an application of the formalism presented in the previous section, we will now 
test the standard cold dark matter (CDM) model of structure formation for consistency with the SP91 
CMB experiment and the Lauer- Postman bulk flow experiment. 

Let El be the Lauer- Postman (LP for short) measurement of bulk flows of galaxies in a 150/i~^Mpc 
sphere (Lauer and Postman, 1993). Let E2 be the 1991 South Pole CMB anisotropy experiment, SP91 
for short (Gaier et al. 1992). Let E3 be the COBE DMR experiment (Smoot et al. 1992). All of these 
experiments probe scales that are well described by linear perturbation theory, and so as long as the 
initial fluctuation are Gaussian, the expected results of the experiments can be expressed simply as 
integrals over the power spectrum of the matter perturbation: 




W^{k)P{k)dk. 



Here Ssp and Scobe are the mean-square temperature fluctuations measured by the experiments, and 
Sip = {v/cf' is the squared bulk flow. The corresponding window functions Wip, Wgp and Wcobe are 
derived in Appendix A, and plotted in Figure 2. These window functions assume that the initial 
perturbations were adiabatic, that f7 = 1, and that recombination happened in the standard way, i.e. 
a last-scattering surface at z w 1000. The SP91 window function is to be interpreted as a lower limit 
to the true window function, as it includes contributions only from the Sachs- Wolfe effect, not from 
Dopplcr motions or intrinsic density fluctuations of the surface of last scattering. 

Now let us turn to the probability distributions for the random variables sip, Ssp and Sc- The 
standard CDM model with power-law initial fluctuations oc fc" predicts a power spectrum that is well 
fitted by (Bond and Efstathiou 1984) 

(l + N+ (69)1-5 + (cg)2]^-i3) 

where a = 6.4, b = 3.0, c = 1.7 and q = {lh~^Mpc)k/r. For the simplest model, T = h, but certain 

additional complications such as a non-zero cosmological constant A and a non-zero fraction Cli, of hot 
dark matter can be fitted with reasonable accuracy by other values of F (Efstathiou, Bond and White 
1992). Thus the model has three free parameters: n, F and the overall normalization A. Integrating 
the power spectrum against the three window functions yields the values of given in Table 2. The 
two rightmost columns contain the quotients cip/cgp and cip/cc, respectively. 

Table 2. Expected r.m.s. signals for CDM power spectrum with A = (lft,~-'^Mpc)^: 



n F LP SP91 COBE LP/SP91 LP/COBE 



1 


0.5 


9.2 X 10- 


-7 


1.6 X 10- 


-8 


2.0 X 10" 


-8 


56.7 


45.1 


0.7 


0.5 


1.7 X 10- 


-6 


2.9 X 10- 


-8 


5.1 X 10" 


-8 


57.2 


32.7 


2 


0.5 


1.4 X 10- 


-7 


2.6 X 10- 


-9 


1.2 X 10" 


-9 


53.9 


112.2 


1 


0.1 


1.4 X 10" 


-6 


2.3 X 10" 


-8 


4.3 X 10" 


-8 


57.9 


31.9 


1 


10 


2.3 X 10" 


-7 


4.1 X 10" 


-9 


4.6 X 10" 


-9 


55.9 


49.6 
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As can be seen, the dependence on F is quite weak, and the quotient cip/cgp is quite insensitive to the 
spectral index n as well. Let us for definiteness assume the canonical values n = 1 and F = 0.5 in 
what follows. 

These values c, would be the average values of the probability distributions for Sgp and sip if 
there where no experimental noise. We will now model the full probability distributions of the three 
experiments, including the contribution from experimental noise. 

For a bulk flow experiment, the three components v^, Vy and of the velocity vector v are 
expected to be independent Gaussian random variables with zero mean, and 

(|v|')=c,p. (8) 

However, this is not quite the random variable sip that we measure, because of errors in distance 

estimation, etc. Denoting the difference between the observed and true bulk velocity vectors by e, let 
us assume that the three components of e are identically distributed and independent Gaussian random 
variables. This should be a good approximation, since even if the errors for individual galaxies are not, 
the errors in the average velocity e will be approximately Gaussian by the Central Limit Theorem. 
Thus the velocity vector that we measure, v + e, is also Gaussian, being the sum of two Gaussians. 
The variable that we actually measure is sip = |v + ep, so 

sip = ^{cip + Vip)xl, (9) 

where xi has a chi-squared distribution with three degrees of freedom, and Vip is the variance due to 

experimental noise, i.e. the average variance that would be detected even if the true power spectrum 
were P{k) = 0. The fact that the expectation value of the detected signal Sip (which is usually referred 
to as the uncorrected signal in the literature) exceeds the true signal cip is usually referred to as error 
bias (LP; Strauss Cen & Ostriker 1993 hereafter SCO). Error bias is ubiquitous to all experiments of 
the type discussed in this paper, including CMB experiments, since the measured quantity is positive 
definite and the noise errors contribute squared. In the literature, experimentally detected signals are 
usually quoted after error bias has been corrected for, i.e. after the noise has been subtracted from 
the uncorrected signal in For LP, the uncorrected signal is 807 km/s, whereas the signal quoted after 
error bias correction is 689 km/s. 

For the special case of the LP experiment, detailed probability distributions have been computed 
using Monte-Carlo simulations (LP, SCO), which incorporate such experiment-specific complications 
as sampling errors, asymmetry in the error ellipsoid, etc. To be used here, such simulations would 
need to be carried out for each value of cip under consideration. Since the purpose of this section is 
merely to give an example of the test formalism, the above-mentioned x^-approximation will be quite 
sufficient for our needs. 

For the SP91 nine-point scan, the nine true values ATi/T are expected to be Gaussian random 
variables that to a good approximation are independent. They have zero mean, and 

(|Ar,/r|^) = c,^. (10) 

Denoting the difference between the actual and observed values by Si, we make the standard assumption 
that these nine quantities are identically distributed and independent Gaussian random variables. Thus 
the temperature fluctuation that we measure at each point, ATi/T + 6, is again Gaussian, being the 
sum of two Gaussians. The variable that we actually measure is 

1 ^ / AT- \ ^ 1 
^^p= Ql2y~jr + ^i) = 9 (cap + ^«p)X9> (11) 

where xi has a chi-squared distribution with nine degrees of freedom, and Vsp is the variance due to 
experimental noise, the error bias, i.e. the average variance that would be detected even if the true 
power spectrum were P{k) = 0. 
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Wc will use only the signal from highest of the four frequency channels, which is the one likely 
to be the least affected by galactic contamination. Again, although Monte-Carlo simulations would 
be needed to obtain the exact probability distributions, we will use the simple x^-approximation here. 
In this case, the main experiment-specific complication is the reported gradient removal, which is a 
non-linear operation and thus does not simply lead to a x^-distribution with fewer degrees of freedom. 

The amplitude of the COBE signal can be characterized by the variance in AT/T on an angular 
scale of 10°. This number can be estimated from the COBE data set as Sc = (jfoo = ((H-O ± 1.8) x 
10~^)^ (Smoot et al. 1992). The uncertainty in this quantity is purely due to instrument noise, and 
contains no allowance for cosmic variance. We must fold in the contribution due to cosmic variance 
in order to determine the probability distribution for Sc- We determined this probability distribution 
by performing Monte-Carlo simulations of the COBE experiment. We made simulated COBE maps 
with a variety of power spectra (including power laws with indices ranging from to 3, as well as 
delta- function power spectra of the sort described in Section 4) . Wc included instrumental noise in the 
maps, and excluded all points within 20° of the Galactic plane. By estimating Sc from each map, we 
were able to construct a probability distribution corresponding to each power spectrum. In all cases, 
the first three moments of the distribution were well approximated by 

Ml = (Sc) = Cc, 

112 = (sl) - {scf < 0.063c^ + 1.44 X 10-^\ (12) 
H3 = {sD = O.OOQc^. 

Furthermore, in all cases the probability distributions were well modeled by chi-squared distributions 
with the number of degrees of freedom, mean, and offset chosen to reproduce these three moments. 
Note that the magnitude of the cosmic variance depends on the shape of the power spectrum as well 
as its amplitude. The inequality in the above expression for /i2 represents the largest cosmic variance 
of any of the power spectra we tested. Since we wish to set conservative limits on models, we will 
henceforth assume that the cosmic variance is given by this worst-case value. Thus we arc assuming 
that the random variable (sc — so)/As has a chi-squared distribution with S degrees of freedom, where 

So = Ml - 2M2/M3, 

As = /I3/4M2, (13) 

S = 8^il/^il 

The results obtained using these three probability distributions are summarized in Tables 3abc. 
In 3a and 3b, N = 2, and the question asked is whether LP is consistent with COBE and SP91, 
respectively. In Table 3c, TV = 3, and we test all three experiments for consistency simultaneously. 
In each case, the optimum normalization (proportional to the entries labeled "Signal") is different, 
chosen such that the consistency probability for the experiments under consideration is maximized. 
As can be seen, the last two tests rule out CDM at a significance level of 95%, i.e. predict that in an 
ensemble of universes, results as extreme as those we observe would be obtained less than 5% of the 
time. Note that using both COBE and SP91 to constrain LP yields a rejection that is no stronger than 
that obtained when ignoring COBE. In the latter case, the best fit is indeed that with no cosmological 
power at all, which agrees well with the observation of SCO that sampling variance would lead LP to 
detect a sizable bulk flow (before correcting for error bias) even if there where none. 

Table 3a. Are LP and COBE consistent with CDM? 





LP 


COBE 


Combined 


Noise 


420km/s 


9.8 /iK 




Signal 


169 km/s 


33.8 /iK 




Noise-|- Signal 


453km/s 


35.2 




Detected 


807 km/s 


35.2 




V 


0.046 


1.00 


0.046 


P{v < v) 


0.046 


1.00 


0.19 
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Table 3b. Are LP and SP91 consistent with CDM? 





LP 


SP91 


Combined 


Noise 


420km/s 


26.4 /zK 




Signal 


km/s 


OfjK 




Noise+ Signal 


420kni/s 


26.4 /xK 




Detected 


807 km/s 


19.9 AtK 




77 


0.023 


0.35 


0.0079 


P{t] < fi) 


0.023 


0.35 


0.046 



Table 3c. Are LP, SP91 and COBE all consistent with CDM? 





LP 


SP91 


COBE 


Combined 


Noise 


420 km/s 


26.4 


9.8 mK 




Signal 


168 km/s 


26.9 /xK 


33.8 /iK 




Noise+Signal 


452 km/s 


37.7 nK 


35.1 /xK 




Detected 


807 km/s 


19.9 AiK 


35.2 yuK 




V 


0.046 


0.039 


0.97 


0.0017 


P{t] < fi) 


0.046 


0.039 


0.97 


0.046 



Tables Sabc show the consistency probability calculations. The first line in each table gives the 
experimental noise, i. e. the detection that would be expected in the absence of any cosmological signal. 
The second line is the best-fit value for the cosmological signal c, the value that maximizes the combined 
consistency probability in the lower right corner of the table. The third line contains the expected value 
of an experimental detection, and is the sum in quadrature of the two preceding lines. The fourth line 
gives the goodness- of- fit parameter for each of the experiments, i.e. the probability that they would 
yield results at least as extreme as they did. The rightmost number is the combined goodness-of-fit 
parameter, which is the product of the others. The last line contains the consistency probabilities, the 
probabilities of obtaining goodness-of-fit parameters at least as low as those on the preceding line. 



4. ALLOWING ARBITRARY POWER SPECTRA 

In this section, we will derive the mathematical formalism for testing results from multiple ex- 
periments for consistency, without making any assumptions whatsoever about the power spectrum. 
This approach was pioneered by Juszkiewicz et al. (1987) for the case N = 2. Here we generalize the 
results to the case of arbitrary N. Despite the fact that the original optimization problem is infinite- 
dimensional, the necessary calculations will be seen to be of a numerically straightforward type, the 
case of independent constraints leading to nothing more involved than numerically solving a system 
of n coupled non-linear equations. After showing this, we will discuss some inequalities that provide 
both a good approximation of the exact results and a useful qualitative understanding of them. 
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4.1. The Optimization Problem 



Let us consider N = n + 1 experiments numbered 0, 1, n that probe the cosmological power 
spectrum P{k). We will think of each experiment as measuring some weighted average of the power 
spectrum, and characterize an experiment Ei by its window function Wi{k) as before. 

Purely hypothetically, suppose we that we had repeated the same experiments in many different 
locations in the universe, and for all practical purposes knew the quantities Ci,...,c„ exactly. Then 
for which power spectrum P{k) would cq be maximized, and what would this maximum be? If we 
experimentally determined cq to be larger than this maximum value, our results would be inconsistent, 
and we would be forced to conclude that something was fundamentally wrong either with our theory 
or with one of the experiments. In this section, we will solve this hypothetical problem. After this, it 
will be seen that the real problem, including cosmic variance and experimental noise, can be solved in 
almost exactly the same way. 

The extremal power spectrum we are looking for is the solution to the following linear variational 
problem: 

Maximize 

/ P{k)Wo{k)dk (14) 
Jo 

subject to the constraints that 

j;^ P{k)Wi{k)dk = a for = l,...,n, 
P{k) > for all A; > 0. 



This is the infinite-dimensional analogue of the so called linear programming problem, and its 
solution is quite analogous to the finite-dimensional case. In geometrical terms, we think of each power 
spectrum as a point in the infinite dimensional vector space of power spectra (tempered distributions 
on the positive real line, to be precise), and limit ourselves to the subset O, of points where all the 
above constraints are satisfied. We have a linear function on this space, and we seek the point within 
the subset Q where this function is maximized. We know that a differentiable functional on a bounded 
region takes its maximum either at an interior point, at which its gradient will vanish, or at a boundary 
point. In linear optimization problems like the one above, the gradient (here the variation with respect 
to P, which is simply the function Wq) is simply a constant, and will never vanish. Thus any maximum 
will always be attained at a boundary point. Moreover, from the theory of linear programming, we 
know that if there are n linear constraint equations, then the optimum point will be a point where 
all but at most n of the coordinates are zero. It is straightforward to generalize this result to our 
infinite-dimensional case, where each fixed k specifies a "coordinate" P{k), and the result is that the 
solution to the variational problem is of the form 

n 

P{k) = Y,P^Sik-ki). (16) 

This reduces the optimization problem from an infinite-dimensional one to a 2n-dimensional one, where 
only the constants pi and ki remain to be determined: 
Maximize 



subject to the constraints that 



Y,PjWo{kj) (17) 



jyj=iPjWi(kj) = a for i = 1, ...,n, ^^^^ 
Pi>0 for i = 1, n. 
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This problem is readily solved using the method of Lagrange multipliers: defining the Lagrangian 



^PjWi{kj) - Ci 



(19) 



and requiring that all derivatives vanish leaves the following set of 3n equations to determine the 3n 
unknowns ft, ki and A^: 

r n 

Wo{ki)-Y,>^jWj{ki) = 0, 



w;,ih)-J2\,w;ih) 



Pi = 0, 



(20) 



Ci-J2PjWi(''j)=^- 

Introducing matrix notation by defining the fc^-dependent quantities Aij = Wj{ki), Bij = Wj(fci), 
flj = WQ{ki) and hi = H^o(A;i) brings out the structure of these equations more clearly: If 0, then 



■ AX 


= a, 


BX 


= b, 




= c. 



(21) 



If A and B are invertible, then eliminating A from the first two equations yields the following system 
of n equations to be solved for the n unknowns fci, A;„: 



A-^a = B-^h. 



(22) 



Although this system is typically coupled and non-linear and out of reach of analytical solutions for 
realistic window functions, solving it numerically is quite straightforward. A useful feature is that once 
this system is solved, a, b, A and B are mere constants, and the other unknowns are simply given by 
matrix inversion: 



X = A- 



(23) 



Since the non-linear system (22) may have more than one solution, all solutions should be substituted 
back into (17) to determine which one is the global minimum. Furthermore, to make statements about 
the solution to our original optimization problem (14), we need to consider also the case where one 
or more of the n variables pi,...,Pn vanish. If exactly m of them are non- vanishing, then without loss 
of generality, we may assume that these are the first m of the n variables. Thus we need to solve the 
maximization problem (17) separately for the cases where P{k) is composed of n delta functions, n — 1 
delta functions, etc., all the way down to the case where P{k) is single delta function. These solutions 
should then be substituted back into (14) to determine which is the global minimum sought in our 
original problem. Thus the solutions depend on the window functions Wj and the signals q in the 
following way: 

•Prom the window functions alone, we can determine a discrete and usually finite number of 

candidate wavenumbers k where delta functions can be placed. 
•The actual signals Cj enter only in determining the coefiicients of the delta functions in the 

sum, i.e. in determining what amount of power should be hidden at the various candidate 

wavenumbers. 
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If wc have found an the optimal sohition, then a small change in the signal vec;tor c will typically 
result in a small change in p and no change at all in the number of delta functions in P{k) or their 
location. If c is changed by a large enough amount, the delta functions may suddenly jump and/or 
change in number as a different solution of (17) takes over as global optimum or one of the coefficients 
Pi becomes negative, the latter causing the local optimum to be rejected for constraint violation. Thus 
within certain limits, we get the extremely simple result that for the optimal power spectrum P{k), 

POO 

Co = / P{k)WQ{k)dk = (A-^a) ■ c. (24) 
Jo 

Thus within these limits, cq depends linearly on the observed signal strengths Cj. This is exactly 
analogous to what happens in linear programming problems. 

4.2. A Useful Inequality 

Before proceeding further, we will attempt to provide a more intuitive understanding of the results 
of the previous section, and show how to determine how complicated a calculation is justified. For the 
special case of only a single constraint, i.e. n = 1, we obtain simply P{k) = pi6{k — ki), where ki is 
given by 

W^{ki)Wi{ki) = Wi{ki)Wo{ki). (25) 
For the case of n constraints, let us define the functions 

Wiikf' 

Then we see that for n = 1, ki is simply the wavenumbcr for which the function /i is maximized, and 
that the maximum signal possible is simply cq = fi{ki). Thus the maximum signal in experiment 
that is consistent with the constraint from experiment i is obtained when the power is concentrated 
where the function fi is large. In other words, if we want to explain a high signal cq in the face 
of low signals in several constraining experiments, then the best place to hide the necessary power 
from the i*^ experiment is where /, takes its maximum. These functions are plotted in Figure 4 for 
the experiments discussed in Section 3, the optimization problem being the search for the maximum 
LP signal that is consistent with the constraints from SP91 and COBE. For illustrative purposes, we 

1 /2 

here assume that Csp and Cg are known exactly, and given by the detected signals Ssp ~ 19.9 /xK and 

1/2 

Sc ~ 33.8 /iK (we will give a proper treatment of cosmic variance and noise in the following section). 

Using the n = 1 constraint for each constraining experiment separately, the smallest of the functions 
thus sets an upper limit to the allowed signal cq = Csp. Thus the limit is given by the highest point in 
the hatched region in Figure 4, i.e. 

Co < c^L = sup min/,(fc). (26) 
k » 

We sec that using the SP91 constraint alone, the LP signal would be maximized if all power were at 
k « (940Mpc)~^. Since this flagrantly violates the COBE constraint, the best place to hide the power 
is instead at fc « (100Mpc)~^ 

By using the above formalism to impose all the constraints at once, the allowed signal obviously 
becomes lower. If the constraints are equalities rather than inequalities, then this stronger limit can 
never lie below value at (fc* ~ 250Mpc)~^, where /sp(fc*) = fc{k*), since this is the signal that 
would result from a power spectrum of the form P{k) cx S{k — fc»). Thus for the particular window 
functions in our example, where the constraint from the n = 2 calculation cannot be more than a 
factor /sp(80Mpc)//sp(250Mpc) « 1.05 stronger than the simple n = 1 limits, the latter are so close 
to the true optimum that they arc quite sufficient for our purposes. If the constraints are upper limits 
rather than equalities, then the limit on cq is more relaxed, and is always the uppermost point in the 
hatched region, i.e. Cmax- 
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4.3. Including Noise and Cosmic Variance 

To correctly handle cosmic variance and instrumental noise, we need to use the formalism devel- 
oped in Section 2. Thus given the probability distributions for the various experimental results Sj, we 
wish to find the power spectrum for which the consistency probability rj is maximized. This optimiza- 
tion problem, in which all experiments are treated on an equal footing, will be seen to lead directly to 
the asymmetric case above where the signal in one is maximized given constraints from the others. For 
definitcncss, wc will continue using the example with the LP, SP91 and COBE experiments. As seen 
in Section 2, the source of the low consistency probabilities is that sip is quite high when compared 
to Ssp and Sc- Thus it is fairly obvious that for the power spectrum that maximizes the consistency 
probability, we will have sip > (sip), whereas Sgp < (ssp) and Sc < (sc), so we can neglect power spectra 
that do not have this property. Let us first restrict ourselves to the subset of these power spectra for 
which cip = D and Cc = E, where D and E are some constants. Then these power spectra all predict 
the same probability distributions for sip and Sc- The consistency probability rj is clearly maximized 
by the power spectrum that maximizes (sgp) , and this will be a linear combination of one or two delta 
functions as shown in Section 4.1. The key point is that since the locations of these delta functions 
are independent of D and E (within the range discussed in 4.1), the infinite-dimensional optimization 
problem reduces to the following two simple steps: 

1. Solve for the optimal number of delta functions m and their locations ki as described in section 

4.1 

2. Find the m coefficients pi for which the power spectrum P{k) = J2iLiPiS{k — ki) maximizes 

the consistency probability. 

4.4. Power spectrum independent constraints on LP, SP91 and COBE 

When applying the above consistency test to the LP, SP91 and COBE experiments, we obtain 
exactly the same consistency probability as in Table 3b. The reason for this is that the optimal 
normalization turns out to be zero. This will obviously change if the LP error bars become smaller in 
the future. Thus dropping the CDM assumption does not improve the situation at all, which indicates 
that main source of the inconsistency must be something other than the CDM model. 

In anticipation of future developments, consistency probabilities where also computed for a number 
of cases with less noise in the LP experiment. Comparing only LP and SP91, the optimum power 
spectrum has a delta function at k Ri (941Mpc)~^. When including all three experiments, treating 
the COBE and SP91 constraints as upper limits, the optimum power spectrum has a single delta 
function at fc « (79Mpc)^^, so the addition of COBE strengthens the constraint only slightly, due to 
the flatness of fsp in Figure 4. Interestingly, for all these cases with smaller LP error bars, consistency 
probabilities were found to be almost as low when allowing arbitrary power spectra as for the CDM 
case. This is again attributable to the flatness of fgp, since weighted averages of a flat function are 
fairly independent of the shape of the weight function (here the power spectrum). 
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5. DISCUSSION 

We have developed a formalism for testing multiple cosmological experiments for consistency. As an 
example of an application, we have used it to place constraints on bulk flows of galaxies using the 
COBE and SP91 measurements of fluctuations in the cosmic microwave background. It was found 
that taken at face value, the recent detection by Lauer and Postman of a bulk flow of 689 km/s on 
scales of 150/i~^Mpc is inconsistent with SP91 within the framework of a CDM model, at a significance 
level of about 95%. However, interestingly, this cannot be due solely to the CDM assumption, since 
the LP result was shown to be inconsistent with COBE and SP91 at the same significance level even 
when no assumptions whatsoever were made about the power spectrum. This leaves four possibilities: 

1. The window functions are not accurate. 

2. Something is wrong with the quoted signals or error bars for at least one of the experiments. 

3. The observed fluctuations cannot be explained within the framework of gravitational instability 

and the Sachs- Wolfe effect. 

4. The random fields are not Gaussian, 

Case 1 could be attributed to a number of effects: HQ ^ 1, then both the calculation of the Sachs- Wolfe 
effect (which determines Wgp and Wcobe) and the growth of velocity perturbations (which determines 
Wip) arc altered. If the universe became reionized early enough to rescatter a significant fraction of 
all CMB photons, then small scale CMB anisotropics were suppressed, which would lower Wgp. A 
quantitative treatment of these two cases will be given in a future paper. Other possible causes of 1 
include a significant fraction of the density perturbations being isocurvature (entropy) perturbations 
or tensor- mode perturbations (gravity waves). Apart from these uncertainties, we have made several 
simplifying assumptions about the window functions for LP and SP91. To obtain more accurate 
consistency probabilities than those derived in the present paper, a more accurate LP window function 
should be used that incorporates the discreteness and the asymmetry of the sample of Abell clusters 
used. This can either be done analytically (Fcldman & Watkins 1993) or circumvented altogether by 
performing Monte-Carlo simulations like those of LP or SCO, but for the whole family of power spectra 
under consideration. 

As to case 2, there has been considerable debate about both the LP and the SP91 experiments. A 
recent Monte-Carlo Simulation of LP by SCO basically confirms the large error bars quoted by LP. As 
is evident from the flatness of LP curve in Figure 3, it will be impossible to make very strong statements 
about inconsistency until future experiments produce smaller error bars. With the SP91 experiment, 
a source of concern is the validity of using only the highest of the four frequency channels to place 
limits, even though it is fairly clear that the other three channels suffer from problems with galactic 
contamination. The situation is made more disturbing by the fact that a measurement by the balloon- 
borne MAX experiment (Gundersen et al. 1993) has produced detections of degree-scale fluctuations 
that that are higher than those seen by SP91, and also higher than another MAX measurement 
(Meinhold et al. 1993). On the other hand. SP91 has been used only as an upper limit in our 
treatment, by including only the Sachs- Wolfe effcc;t and neglecting both Dopplcr contributions from 
peculiar motions of the surface of last scattering and intrinsic density fluctuations at the recombination 
epoch. If these effects (which unfortunately depend strongly on parameters such as h and fif,) where 
included, the resulting constraints would be stronger. 

Case 3 might be expected if the universe underwent a late-time phase transition, since this could 
generate new large-scale fluctuations in an entirely non-gravitational manner. 

In the light of the many caveats in categories 1 and 2, the apparent inconsistency between LP 
and SP91 (Jaffe et al. 1993) is hardly a source of major concern at the present time, and it does not 
appear necessary to invoke 3 or 4. However, we expect the testing formalism developed in this paper 
to be able to provide many useful constraints in the future, as more experimental data is accumulated 
and error bars become smaller. 

The authors wish to thank Michael Strauss, Bhuvnesh Jain, Douglas Scott, Martin White, and 
Joseph Silk for many useful comments. 
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APPENDIX: WINDOW FUNCTIONS 

The results of CMB anisotropy experiments can be conveniently described by expanding the 
temperature fluctuation in spherical harmonics: 

AT °° ' 

^(f) = E E (27) 

1=2 m=-l 

(The monopole and dipole anisotropics have been removed from the above expression, since they are 
unmeasurablc.) If the fluctuations are Gaussian, then each coefficient aim is an independent Gaus- 
sian random variable with zero mean (Bond and Efstathiou 1987). The statistical properties of the 
fluctuations are then completely specified by the variances of these quantities 

Ci = {\ai^f). (28) 

(The fact that the variances are independent of m is an immediate consequence of spherical symmetry.) 
Different CMB experiments are sensitive to different linear combinations of the Ci's: 

oo 

S = J2FiCi, (29) 

1=2 

where S is the ensemble-averaged mean-square signal in a particular experiment, and the "filter func- 
tion" Fi specifies the sensitivity of the experiment on different angular scales. The filter functions for 
COBE and SP91 are 

(co6e) _ (^[±11 _-<7^(;-Fi) = 
' ~ 47r 

; (30) 

m.=—l 

where Hq is a Struve function. CTc = 4.25° and (Tg = 0.70° are the r.m.s. beamwidths for the two 
experiments, and a = 1.5° is the amplitude of the beam chop (Bond et al. 1991, Dodelson and Jubas 
1993, White et al. 1992). 

For Sachs- Wolfe fluctuations in a spatially flat Universe with the standard ionization history, the 
angular power spectrum C; is related to the power spectrum of the matter fluctuations in the following 
way (Peebles 1984, Bond and Efstathiou 1987): 



Ci = — dkP{k)jt{k). (31) 
t^Tq Jo 

Here tq is the conformal time at the present epoch, and 

jl{k)= J jl{kT)V{kT)dT, (32) 

where ji is a spherical Bcssel function. The visibility function V is the probability distribution for the 
conformal time at which a random CMB photon was last scattered. ji{k) is therefore the average of 
jiikr) over the last scattering surface. We have used the V of Padmanabhan (1993). 

We can combine equations (29), (30), and (31) to get the window functions for the two experiments: 

(_o 

1 I (33) 

1=2 m=-l 
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The mean-square bulk flow inside of a sphere of radius a is (see, e.g., Kolb and Turner 1990) 



However, we must make two corrections to this result before applying it to the LP data. This formula 
applies to a measurement of the bulk flow within a sphere with an inflnitcly sharp boundary. In reality, 
errors in measuring distances cause the boundary of the spherical region to be somewhat fuzzy. If 
we assume that distance measurements are subject to a fractional error e, then the window function 
must be multiplied by er'^'^^'^^ . We have taken e = 0.16, the average value quoted by LP. It shoifld be 
noted that this value varies from galaxy to galaxy in the LP sample, due to the distance estimation 
technique used, and that a more accurate window function that reflects the discrete locations of the 
Abell clusters used in the survey should take this into account. 

The second correction has to do with the behavior of the window function at small k. Equation 
(34) applies to the velocity relative to the rest frame of the Universe. The velocity measured by LP 
is with respect to the CMB rest frame. If there is an intrinsic CMB dipolc anisotropy, then these 
two reference frames differ. Therefore, we must include in equation (29), a term corresponding to the 
intrinsic CMB dipole. This correction was first noticed by Gorski (1991). After applying both of these 
corrections, the LP window function is 




(34) 




2 



(35) 
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k [h-iMpc] 




True LP power [km/s] 



